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FOREWORD 


This  report  summarizes  an  Investigation  of  three  Newton-llke  methods  for  the 
numerical  solution  of  the  small  disturbance  equation  of  non-lifting  transonic 
flow  past  a  parabolic  airfoil.  This  work  was  performed  by  the  authors,  who  are 
members  of  the  Applied  Mathematics  Branch  of  NSWC.  The  authors  wish  to  thank 
Dr.  Gregory  R.  Shubin  of  NSWC  for  his  interest  in  the  work  reported  herein, 
thanking  him  for  his  many  comments  and  suggestions. 

fcsjiWv*' 
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1.  INTRODUCTION 

In  this  paper  we  investigate  the  effectiveness  of  Newton’s  method  and  two 
variants  of  Newton's  method  for  solving  the  nonlinear  difference  equations  that 
approximate  the  transonic  small  disturbance  equation.  In  particular,  we  consider 
the  case  of  non-lifting  flow  over  a  parabolic  airfoil.  Numerical  methods  for 
solving  both  the  small  disturbance  and  the  full  potential  equations  have  been  the 
object  of  a  good  deal  of  research  in  the  last  few  years.  The  references  [1]  - 
[10]  give  an  idea  of  the  range  of  different  methods  that  have  appeared  in  the 
literature. 

Although  Newton's  method  converges  quadratlcally,  requiring  fewer  iterations 
than  linear ly-convergent  methods,  its  effectiveness  for  nonlinear  problems  with 
many  unknowns  is  open  to  question  for  several  reasons.  First,  at  each  iteration 
one  must  typically  solve  a  large  system  of  linear  equations.  The  relative 
success  of  Newton's  method  depends  in  part  on  how  effectively  one  can  solve  the 
linear  system  either  directly  or  iteratively.  A  direct  method  may  require  an 
excessive  amount  of  computer  memory,  and  require  out-of-core  solvers,  and  for 


1.  Bailey,  F.  R. ,  "On  the  Computation  of  Two-  and  Three-dimensional  Steady 
Transonic  Flows  by  Relaxation  Methods,"  Lecture  Notes  in  Physics,  v.  41: 
Progress  in  Numerical  Fluid  Dynamics  (H.  J.  Wirz,  ed.),  Springer-Verlag. 
Berlin,  pp.  1-77,  1975. 

2.  Enquist,  B.  and  Osher,  S.,  "Stable  and  Entropy  Satisfying  Approximations  for 
Transonic  Flow  Calculations,"  Math.  Comp. ,  v.  34,  pp.  45-75,  January  1980. 

3.  Enquist,  B.  and  Osher,  S. ,  "One-sided  Differences  and  Transonic  Flow," 
to  appear,  Proc.  Nat.  Acad.  Sci. 

4.  Jameson,  A.,  "Iterative  Solution  of  Transonic  Flows  Over  Airfoils  and  Wings, 
Including  Flows  at  Mach  1,"  Comm.  Pure  Appl.  Math.,  v.  27,  pp.  283-309,  1974. 

5.  Jameson,  A.,  "Numerical  Solution  of  Nonlinear  Partial  Differential  Equations 
of  Mixed  Type,"  in  Numerical  Solution  of  Partial  Differential  Equations  III 
(B.  Hubbard,  ed.),  SYNSPADE  1975,  pp.  275-320,  Academic  Press,  New  York,  1976. 

6.  Hafez,  M.  M.  and  Cheng,  H.  K. ,  "Convergence  Acceleration  of  Relaxation 
Solutions  for  Transonic  Flow  Computations,"  AIAA  Journal,  v.  15,  no.  3, 
pp.  329-336,  1977. 

7.  Hafez,  M.  M. ,  South,  J. ,  and  Murman,  E. ,  "Artificial  Compressibility  Methods 
for  Numerical  Solutions  of  Transonic  Full  Potential  Equation,"  AIAA  Journal, 
v.  17,  no.  8,  pp.  838-844,  1978. 

8.  Majda,  A.  and  Osher,  S.,  "Numerical  Viscosity  and  the  Entropy  Condition," 

Comm.  Pure  Appl.  Math,  (to  appear) . 

9.  Martin,  E.  D. ,  "A  Fast  Semi-direct  Method  for  Computing  Transonic  Aerodynamic 
Flows,"  Proc.  Second  AIAA  Conf.  Comput.  Fluid  Dyn. ,  Hartford,  pp.  162-174, 

June  1975. 

10.  South,  J.  C.  and  Brandt,  A.,  "Application  of  a  Multi-level  Grid  Method  to 
Transonic  Flow  Calculations,"  ICASE  Report  # 76-8,  Hampton,  Va. ,  1976. 
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problems  of  the  type  we  are  considering  an  iterative  method  such  as  SOR  may  fail 
to  converge.  Secondly,  the  domain  of  convergence  for  Newton's  method  may  be 
relatively  small,  requiring  a  good  initial  guess  for  convergence. 

The  two  variants  of  Newton's  method  we  consider  seem  to  be  more  robust  than 
Newton's  method  with  respect  to  convergence  and  require  less  storage;  moreover, 
one  requires  less  computer  time  to  converge. 

In  section  2  we  discuss  the  formulation  of  the  transonic  small  disturbance 
equation  and  the  Murman-Cole  difference  scheme  for  its  approximation.  In 
section  3  we  discuss  the  two  variants  of  Newton's  method  for  the  nonlinear 
difference  equations.  One  variant  is  a  line-Newton  successive  over-relaxation 
(LINSOR)  method  and  the  other  uses  only  the  lower  triangular  portion  of  the 
Jacobian  with  over-relaxation  of  the  iterates  (LONSOR) . 

The  numerical  performance  of  these  methods  and  Newton's  method  are  compared 
in  section  4.  Briefly,  these  numerical  results  indicate,  for  example,  that  with 
the  proper  choice  of  an  over-relaxation  parameter,  the  LINSOR  method,  when 
compared  to  Newton's  method  is  faster,  needs  a  good  deal  less  storage,  and 
converges  for  a  wider  range  of  starting  guesses.  The  LONSOR  method  also  needs 
less  storage  and  on  a  25  x  25  mesh  satisfies  the  steady  state  error  criterion  in 
a  time  comparable  to  Newton's  method. 

2.  FORMULATION  OF  THE  PROBLEM 

We  consider  the  two  dimensional  small  disturbance  equation 

S($)  =  (K-Cy+D^)^  +  $yy  -  0  (2.1) 

which  models  inviscid  flow  over  a  thin  wing.  Here  d>  is  the  potential  function 
of  the  perturbation  from  the  free  stream  velocity.  The  transonic  similarity 
parameter  K  is  defined  by 

K  -  (l-M2)/62/3 
00 

where  M  is  the  free  stream  Mach  number,  assumed  to  be  close  to  1,  and  6  is  the 
airfoil  thickness  ratio.  For  a  derivation  of  (2.1)  from  the  Euler  equations  see, 
for  example,  [11]. 

BOUNDARY  CONDITIONS.  In  order  to  approximate  boundary  conditions  at 
infinity  the  airfoil  is  enclosed  within  a  sufficiently  large  rectangle  R.  For 

_  i 

11.  Nahamira,  S.,  Computational  Methods  in  Engineering  and  Science,  Wiley- 
Interscience,  New  York,  1977. 
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lifting  flows,  <p  at  the  boundary  3R  of  R  is  defined  in  terms  of  the  circulation. 
For  nonlifting  flows  we  set  <p  «  0  on  }R. 

At  the  boundary  of  the  airfoil  corresponding  to  a  point  (x,  0)  it  is  required 
that  0+)  ■  F^(x)  and  4>y(x,  0  )  ■  FVx)  where  F+(x)  and  F_(x)  describe  the 

upper  and  lower  boundaries  of  the  airfoil  respectively.  We  consider  a  parabolic 
airfoil  where  F+(x)  and  F_(x)  are  given  by  25x(l-x)  and  -25x(l-x)  respectively. 
These  conditions  are  applied  on  the  mean  surface  y  ■  0  and  approximate  the 
condition  that  the  flow  at  the  airfoil  be  tangent  to  the  airfoil. 

FLOW  OF  MIXED  TYPE.  Equation  (2.1)  together  with  the  boundary  conditions 
constitute  a  nonlinear  boundary  value  problem  for  (f>.  Basically  the  type  of  flow 
at  a  given  point  is  regulated  by  the  sign  of  K  -  (y+l)$  .  For  a  given  6  when  M 
is  sufficiently  close  to,  but  less  than,  one,K  -  (y+l)<|>  is  negative  in  certain 
regions  above  and  below  the  airfoil.  In  these  regions  the  flow  is  supersonic 
and  is  bounded  by  a  sonic  line  and  weak  shock.  Outside  this  region  K  -  (y+l)<J>x 
is  positive  and  the  flow  is  subsonic. 

Different  numerical  techniques  have  been  devised  to  approximate  the  solution 
We  now  describe  a  finite  difference  scheme  and  iterative  methods  for  its 
solution. 


FINITE  DIFFERENCE  EQUATIONS.  A  prevalent  method  for  discretizing  (2.1)  is 
the  Murman-Cole  difference  scheme  introduced  in  [12]  and  later  modified  to  more 
accurately  capture  shocks  [13].  A  more  recent  modification  which  excludes 
expansion  shocks  is  due  to  Enquist  and  Osher  [2,  3]. 

In  the  Murman-Cole  scheme  directionality  is  introduced  by  using  centered 


differences  for  <t>x  and  in  the  subsonic  region  and  backward  differences  in  the 
supersonic  region.  This  procedure  is  implemented  as  follows.  We  partition  the 


rectangle  R  with  a  uniform  mesh  x.  -  ih,  0  1  i  <  n+1,  y.  -  jh,  0  t  j  1  m+1.  We 

(k)  th  1  J 

define  to  be  the  approximation  to  $(x^,  y^)  defined  by  some  iterative 

scheme.  At  the  mesh  point  (x^,  y^)  define 


12.  Murman,  E.  M.  and  Cole,  J.  D.,  "Calculation  of  Plane  Steady  Transonic  Flows," 
AIAA  Journal,  v.  9,  no.  1,  pp.  114-121,  1971. 

13.  Murman,  E.  M. ,  "Analysis  of  Embedded  Shock  Waves  Calculated  by  Relaxation 
Methods,"  Proceedings  of  AIAA  Computational  Fluid  Dynamics  Conference, 

Palm  Springs,  Ca. ,  pp.  27-40,  July  1973. 
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'ij 


and 


(♦S'”  -  ♦&,)«*« 


Hy  =  K  -  (Y+l)  ^5_1) 


then  E . .  and  H. .  are  discrete  versions  of  the  term  K  -  (y+1)$  . 
il  ij  x 

is  judged  to  be  elliptic  at  (x^,  y^)  and  is  differenced 

E,,524>. .  +  624> . ,  =  0 
i3  x  i]  y  ij 


where 


and 


5x*ij  =  (*i+l,j  "  2*ij  + 


) / (Ax)' 


62<f>.  . 
y  ij 


(*i,j+i  ~  2*ij  + 


Pi,j-1 


)/(Ay)‘ 


If  E±j  >  0,  (2.1) 


If  Ey  <  0  and  <  0  (2.1)  is  hyperbolic  at  (x^  y^)  and  differenced  as 

.2 , 


hiis**i-i,3  +  Vij 


o. 


If  E^  <  0  and 


>  0  the  point  (x^,  y^)  is  regarded  as  lying  on  a  sonic  line 


and  (2.1)  at  (x^  y^)  is  differenced  as 


62*,, 

y  ij 


o. 


In  the  case  that  E^  >  0  and  <  0  the  point  (x^,  y^)  lies  at  a  shock  and  the 
original  Murman-Cole  formulation  was  to  use  centered  differencing.  Later  this 
was  modified  to  incorporate  a  "shock  point"  operator  given  by 


<E1J  +  V6K*ij  +  fyil 


0. 


In  this  paper  we  use  centered  differencing  at  a  shock  point  since  we  are  not 
so  interested  in  the  accuracy  of  the  shock  strength  or  shock  location  but  rather 
in  the  efficient  solution  of  the  finite  difference  equations.  Finite  difference 
approximations  for  the  specification  of  $  on  the  airfoil  boundary  are  actually 
applied  on  the  rows  of  mesh  points  (x^,  yj)  immediately  above  and  below  the  mean 
surface  of  the  airfoil.  There  the  given  expression  for  is  incorporated  by 


.  In  if1  t'lirtfritk'  Ziii.fr.aai 
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»  Si. 

yy  Ay 


(♦ 


LJ±1 


.) 


-*.i**- : .  *  | 

Ay  y 1 y=0 


above  the  wing  and 


yy 


i[< 


(o 


y'y=0 


Ay 


below  the  wing. 

3.  ITERATIVE  METHODS 


We  write  the  nonlinear  mn  x  tan  system  of  equations  as  S($>)  =  0.  A  general 
iterative  scheme  for  the  solution  of  S(<j>)  =  0  is 


N(ft(ft+1 


A  -  -as (**> 


(3.1) 


where  N(if  )  is  an  mn  x  in  matrix  and  a  is  a  relaxation  parameter.  This  iterative 
method  may  be  viewed  as  a  discretization  of  the  time-dependent  equation 

N(<Kt)H(t)  +  S(<Kt)) 


0. 


(3.2) 


th 


In  the  case  that  N  is  the  Jacobian  matrix  J  whose  (i,j) —  entry  is  given  by 
the  derivative  of  the  i—  equation  with  respect  to  the  j—  unknown,  the  iterative 
scheme  (3.1)  is  a  damped  Newton's  method  for  the  system  (3.2);  in  this  case, 

J  has  the  block  form 


•  •  •  • 

B  ,C  .D  ,E  , 
n-1  n-1  n-1  n-1 


The  n  blocks  C^,  D^,  E^  are  m  x  m  matrices.  The  matrix  is  tridiagonal  and 
the  blocks  C^,,  E^  are  diagonal.  The  blocks  B.^  have  non-zero  entries 
corresponding  to  points  (x^,  y^)  in  the  hyperbolic  region. 

LINE-NEWTON  SOR  (LINSOR) .  Despite  the  rapid  convergence  of  Newton's  method 
for  good  initial  guesses,  it  suffers  from  two  major  drawbacks:  its  storage 
requirement  of  ©(m^n)  entries,  due  to  the  fill-in  which  occurs  when  Gaussian 
elimination  is  used  to  solve  the  system  (3.1)  with  N  *  J,  and  its  operation  count 


9 


NSWC  TR  80-271 


of  0(m3n)  operations  per  iteration.  In  order  to  facilitate  the  solution  of  (3.1) 

by  Gaussian  elimination,  we  let  N  be  the  block  lower  triangular  matrix  J*  which 

is  formed  by  deleting  the  blocks  from  the  matrix  J.  If  the  diagonal  blocks 

are  nonsingular,  we  may  then  solve  (3.1)  as  a  series  of  block  solves:  writing 

yk+1  -k  — k  ,  T  T.T 

<j>  -  4>  =  e  =  (e1*  **en) 


S(^)  =  (f^'-fV, 


we  have  a  series  of  n  tridiagonal  systems 


Dlel  =  -fl 

Vi  =  "fi  "  Ciei-1 


(2  <  i  <  n) . 


— k+1 

In  order  to  accelerate  convergence,  we  overwrite  tp  by  a  convex  combination  of 
the  old  and  new  iterates 

,k+l  ,k+l  .  k 

d>±  +  tui4>±  +  (1-ukH.. 

—k+1 

Here  the  relaxation  parameter  satisfies  >  1  for  those  components  of  $ 
lying  in  the  supersonic  region  and  us^  *  1  for  those  components  in  the  sonic  and 
subsonic  regions.  Note  that  this  method  requires  storage  0(m)  and  uses  0(mn) 
operations  per  iteration,  so  that  an  iteration  of  the  LINSOR  method  is  far 
cheaper  than  an  iteration  of  Newton's  method  and  the  LINSOR  method  has  a  much 
smaller  storage  overhead  than  does  Newton's  method. 

One  may  show  that  the  method  above  is  equivalent  to  the  following.  We 
consider  the  difference  equations  approximating  (1.1)  and  choose  a  column  £  in 

the  finite-difference  mesh.  We  consider  the  values  ,  •  •  •  ,4>0  on  column  £  as 

J6X  x,m 

the  only  unknowns  and  consider  only  those  difference  equations  which  involve  the 
$...  In  these  equations,  we  use  the  most  recent  values  of  <t>  .  and  4>  .  when 

they  are  needed.  One  iteration  of  Newton's  method  is  used  to  approximately 

solve  for  the  unknowns  4>.  .  These  values  are  then  relaxed  where  appropriate  and 

Jcj 

updated.  We  then  move  to  column  £+1  and  repeat  the  process.  Starting  at  column  1 
and  continuing  to  column  n  constitutes  a  sweep.  The  grid  is  repeatedly  swept 
until  steady  state  is  reached.  This  technique  may  be  viewed  as  a  type  of  line 
SOR-Newton  sweep. 
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LOWER  NEWTON  SOR  (LONSOR) .  We  may  also  construct  a  new  iteration  matrix  N 
by  modifying  J  in  the  following  simple  manner:  we  eliminate  from  J  all  entries 
above  the  diagonal,  leaving  a  lower  triangular  matrix  with  nonzero  elements  on 
the  diagonal.  We  may  quickly  solve  the  resulting  system  using  forward  substitu¬ 
tion.  Those  $  which  are  in  the  elliptic  region  of  flow  can  be  relaxed  as  soon 
as  they  are  computed.  The  matrix  N  can  now  be  regarded  as  the  Jacobian  matrix 

corresponding  to  the  equations  (3.1),  where  at  the  point  (x . ,  y.),  tj>  .  and 

i  J  >  3 

<f> .  .  .  are  no  longer  treated  as  unknowns,  but  instead  their  most  recent  values 
are  substituted. 

The  numerical  results  for  Newton's  method,  L1NS0R,  and  LONSOR  are  presented 
in  the  next  section. 

4.  RESULTS 

In  this  section,  we  present  computational  results  for  the  problem  stated  in 

Section  3.  The  computational  region  was  the  rectangle  [-1.5,  1.5]  *  [-1.0,  1.0] 

of  the  x-y  plane,  the  wing  being  a  parabolic  airfoil  on  the  segment  [0,  1]  of  the 

x-axis.  The  physical  parameters  used  were  =  0.9  for  the  free-stream  Mach 

number  and  5  between  0.05  and  0.065  for  the  airfoil  thickness  ratio.  The  problem 

was  run  on  uniform  20*20  and  25*25  meshes  of  unknowns.  The  error  was  measured 

by  the  maximum  modulus  of  the  residual  vector,  and  the  stopping  criterion  was  that 

_4 

the  error  be  less  than  10 

The  three  methods  that  were  successfully  tested  are  the  Line-Newton  SOR 
(LINS0R)  method,  the  Lower-Newton  SOR  (LONSOR)  method,  and  the  full  Newton  method. 
For  the  results  reported  in  this  section  an  initial  guess  of  $  =  0  was  used  for 
all  three  methods. 

The  LINS0R  and  LONSOR  methods  require  an  iteration  parameter  to.  It  was 
decided  to  use  over-relaxation  only  on  the  elliptic  region,  choosing  to  =  1  in  the 
hyperbolic  and  parabolic  regions.  In  Figure  1,  we  present  experimental  evidence 
which  shows  that  choosing  a  =  1.65  in  the  elliptic  region  is  close  to  optimal  for 
the  20*20  mesh  with  the  LINS0R  method;  it  was  similarly  determined  that  choosing 
o o  »  1.85  in  the  elliptic  region  is  close  to  optimal  for  the  20*20  mesh  with  the 
LONSOR  method.  Moreover,  the  dependence  of  to  on  mesh  size  does  not  appear  to  be 
too  critical;  the  same  values  were  used  successfully  on  the  25*25  mesh. 

Figure  1  indicates  that  elliptic-region  over- relaxation  successfully 
accelerates  convergence.  However,  it  is  not  a  priori  clear  at  which  iteration  the 
over- relaxation  should  start.  For  example,  when  over-relaxation  was  started 
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immediately  for  the  case  =  0.95,  the  method  diverged;  experimentally  it  was 
found  that  for  =  0.95,  one  could  successfully  start  over-relaxation  after  about 
twelve  iterations.  The  shift  from  the  non-over-relaxed  to  the  over-relaxed 
scheme  causes  the  error  to  increase  momentarily. 

Although  the  LINSOR  and  L0NS0R  methods  were  relatively  simple  to  implement, 
there  was  a  storage  problem  with  the  full  Newton  method.  This  problem  was  due  to 
the  fill-in  of  the  factorization  of  the  sparse  Jacobian  matrix  as  the  product  of 
lower-  and  upper-triangular  banded  matrices.  To  alleviate  this  problem,  the  Yale 
sparse  matrix  package  (YSMP)  with  uncompressed  storage  option  for  nonsymmetric 
matrices  [14]  was  used.  Despite  the  improvement  in  storage  when  using  YSMP  instead 
of  a  banded  solver,  the  full  Newton/YSMP  program  required  a  prohibitive  amount  of 
storage  as  the  number  of  unknowns  increased.  For  instance,  the  Newton/YSMP  program 
could  barely  fit  the  25*25  problem  onto  a  CDC  6500,  while  LINSOR  was  able  to  handle 
problems  of  size  40*40  or  larger. 

In  Figures  2  and  3,  results  are  given  for  the  20*20  and  25*25  meshes, 
respectively.  In  each  figure,  the  maximum  modulus  of  the  residual  vector  is 

plotted  as  a  function  of  the  CP  time  used,  for  each  of  the  three  methods  tested. 

-4 

One  sees  that  LINSOR  has  the  fastest  convergence  in  each  case,  meeting  the  10 
error  criterion  for  the  20*20  and  25*25  matrices  in  9.229  and  22.06  seconds, 
respectively.  The  full  Newton  method  was  next  fastest,  requiring  15.37  and  31.62 
seconds  for  the  two  meshes.  LONSOR  was  the  slowest  of  the  three,  requiring  18.74 
and  35.01  seconds. 

Not  surprisingly,  only  a  few  full  Newton  iterations  were  needed  to  solve  the 
-4 

problem  to  within  10  (6  and  7  iterations  for  the  20*20  and  25*25  cases, 

respectively).  However,  the  time  required  to  solve  the  linear  system  resulting 
from  each  Newton  iteration,  combined  with  the  space  required  to  store  the 
triangular  factors  of  the  coefficient  matrix,  rendered  the  full  Newton  method 
inferior  to  LINSOR. 

A  number  of  other  methods  were  attempted.  Following  ideas  of  Jameson  [5],  a 
fast  Poisson  solver  was  used  after  every  few  iterations  of  the  LINSOR  method;  this 


14.  Eisenstadt,  S.  C.,  Gursky,  M.  C.,  Schultz,  M.  H. ,  and  Sherman,  A.  H. ,  "Yale 
Sparse  Matrix  Package,"  Yale  University  Computer  Science  Department  Research 
Report  #114,  New  Haven,  1977. 
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was  found  to  be  difficult  to  implement  effectively.  The  SOR  method,  when  applied 
to  the  linear  system  arising  from  Newton's  method,  was  found  to  be  "non-robust"  in 
the  sense  that  it  failed  to  work  for  =  0.95.  (See  [15]  for  discussion  of  using 
the  SOR  method  for  mixed  problems.) 

The  block- SOR  and  conjugate  gradient  methods  were  also  used  in  the  attempt  to 
solve  the  linear  system  arising  from  Newton's  method.  The  former  sometimes  showed 
oscillatory  behavior  or  diverged,  while  the  latter  required  a  prohibitively  large 
number  of  conjugate  gradient  iterations  to  solve  the  linear  system  with  sufficient 
accuracy. 

It  is  interesting  to  note  that  using  different  methods  to  solve  the  Murman- 
Cole  difference  equations  may  yield  different  solutions.  Using  a  relatively  thick 
wing  at  =  0.95,  and  a  10x10  grid  of  unknowns,  LINSOR  found  a  solution  whose 
sonic  line  was  uniformly  one  mesh  point  to  the  right  of  that  found  by  the  full 
Newton  method  (see  Figure  4).  Since  the  sonic  lines  were  different,  the  type- 
dependence  of  the  difference  scheme  implies  that  the  two  final  systems  of  nonlinear 
equations  were  different,  and  had  different  solutions. 

Finally,  we  remark  that  using  LINSOR  and  a  continuation  technique  with  respect 
to  M^,  it  is  possible  to  compute  solutions  up  to  and  past  Mach  1  for  the  small 
disturbance  equation. 


15.  Steger,  J.  L.  and  Lomax,  H. ,  "Transonic  Flow  About  Two  Dimensional  Airfoils 
by  Relaxation  Procedures,"  AIAA  Journal,  v.  10,  pp.  49-54,  1972. 
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FIGURE 


1  CHOICE  OF  OPTIMAL  u  FOR  THE  LINSOR  METHOD  USING  A  20*20  MESH  | 
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FIGURE  4  SONIC  PROFILES  OF  DIFFERENT  SOLUTIONS  FOUND  BY  THE 
LINSOR  AND  NEWTON  METHODS  ON  A  10x10  MESH 
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